Data Preparation
#For forecasting, we chose the latest data
trend_hosp_hb <- trend_hb_daily %>%
filter (hb_name == "Scotland") %>%
filter(date >="2021-06-01") %>%
filter(!(is.na(hospital_admissions))) %>%
select(date, hospital_admissions)
# Convert it into a time series
daily_hosp_hb_zoo <- zoo(trend_hosp_hb$hospital_admissions,
order.by=as.Date(trend_hosp_hb$date, format='%m/%d/%Y'))
# Convert it into a time series
daily_hosp_hb_timeseries <- timeSeries::as.timeSeries(daily_hosp_hb_zoo)
Step 1 : Visualize the time series
original_series<-autoplot(daily_hosp_hb_timeseries, ts.colour = '#5ab4ac')+
xlab("Month") +
ylab("Number of People hospitalised")+
#scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
ggtitle("Trend on Hospitalisation") +
color_theme()
ggplotly(original_series)
Step 2 : Identification of model : (Finding d:)
Identify whether the time series is stationary / non stationary we can use ADF Augmented Dickey-Fuller test
adf_test_hosp <- adf.test(daily_hosp_hb_timeseries)
adf_test_hosp
Augmented Dickey-Fuller Test
data: daily_hosp_hb_timeseries
Dickey-Fuller = -1.0883, Lag order = 4, p-value = 0.9209
alternative hypothesis: stationary
The time series is not stationary since we have a high p-value (p-value must be < 0.05). So we apply difference
first_diff_hosp<- diff(daily_hosp_hb_timeseries)
adf_test1_hosp <- adf.test(na.omit(first_diff_hosp))
Warning in adf.test(na.omit(first_diff_hosp)) :
p-value smaller than printed p-value
adf_test1_hosp
Augmented Dickey-Fuller Test
data: na.omit(first_diff_hosp)
Dickey-Fuller = -7.5072, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary
Create a dataframe to compare
adf_data_hosp <- data.frame(Data = c("Original", "First-Ordered"),
Dickey_Fuller = c(adf_test_hosp$statistic, adf_test1_hosp$statistic),
p_value = c(adf_test_hosp$p.value,adf_test1_hosp$p.value))
adf_data_hosp
Initially the p-value is high which indicates that the Time Series is not stationary. So we apply difference 1 time. After the first difference, the p-value < significance level (0.05) So we can conclude that the difference data are stationary. So difference (d = 1)
Other method:
ndiffs(daily_hosp_hb_timeseries)
Let’s plot the First Order Difference Series
Order of first difference
p<- autoplot(first_diff_hosp, ts.colour = '#5ab4ac') +
xlab("Month") +
ylab("HOSPITALIZATION")+
# scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
ggtitle("First-Order Difference Series") +
color_theme()
ggplotly(p)
Step 3: Estimate the parameters (Finding p and q)
For our model ARIMA (p,d,q), we found d = 1, the next step is to get the values of p and q, the order of AR and MA part. Plot ACF and PACF charts to identify q and p respectively.
par(mfrow=c(2,2))
acf_hosp <- acf1(daily_hosp_hb_timeseries, col=2:7, lwd=4)
pacf_hosp <- acf1(daily_hosp_hb_timeseries, pacf = TRUE, col=2:7, lwd=4)
acf_hosp <- acf1(first_diff_hosp, col=2:7, lwd=4)
pacf_hosp <- acf1(first_diff_hosp, pacf = TRUE, col=2:7, lwd=4)
The ACF and PACF plots of the differenced data show the following patterns:
The ACF is sinusoidal and there is a significant spike at lag 3 in the PACF So the data may follow an AR(3) model
The PACF is sinusoidal and there is a significant spike at lag 2 in the ACF So the data may follow an MA(2) model
So we propose three ARMA models for the differenced data: ARMA(p,q) ARMA(3,2), ARMA(3,0) and ARMA(0,2).
That is, for the original time series, we propose three ARIMA models,ARIMA(p,d,q) ARIMA(3,1,2), ARIMA(3,1,0) and ARMA(0,1,2).
Step 4: Build the ARIMA model
Manual ARIMA:
arima_fit_hosp_1 = Arima(daily_hosp_hb_timeseries, order = c(3,1,2))
arima_fit_hosp_2 = Arima(daily_hosp_hb_timeseries, order = c(3,1,0))
arima_fit_hosp_3 = Arima(daily_hosp_hb_timeseries, order = c(0,1,2))
summary(arima_fit_hosp_1)
summary(arima_fit_hosp_2)
summary(arima_fit_hosp_3)
texreg::screenreg(list(arima_fit_hosp_1, arima_fit_hosp_2, arima_fit_hosp_3),
custom.model.names =c("ARIMA(3,1,2)","ARIMA(3,1,0)","ARIMA(0,1,4)"),
center = TRUE,
table = FALSE)
#Function for Automated ARIMA
auto_arima_fit_hosp <- auto.arima(lag(daily_hosp_hb_timeseries,k =lag_value ),
seasonal=FALSE,
stepwise=FALSE,
approximation=FALSE,
trace = TRUE
)
Automated ARIMA
#Lag is used to best fit the model
auto_arima_fit_hosp <- auto.arima(lag(daily_hosp_hb_timeseries),
seasonal=FALSE,
stepwise=FALSE,
approximation=FALSE,
trace = TRUE
)
auto_arima_fit_hosp
Automated ARIMA confirms that the ARIMA(3,1,2) seems good based on AIC
coef<-lmtest::coeftest(auto_arima_fit_hosp)
coef
All coefficients are significant except ar3.
Model Selection Criteria :
ARIMA models with minimum AIC, RMSE and MAPE criteria were chosen as the best models. Based on Akaike Information Criterion (AIC) above, an ARIMA(3, 1, 2) model seems best.
Step 5: Check for Diagnostics
Let’s plot the diagnostics with the results to make sure the normality and correlation assumptions for the model hold. If the residuals look like white noise, proceed with forecast and prediction, otherwise repeat the model building.
res <-checkresiduals(auto_arima_fit_hosp, theme = color_theme())
res
The ACF plot of the residuals from the ARIMA(3,1,2) model shows that all auto correlations are almost within the threshold limits, with residuals. A portmanteau test (Ljung-Box test) returns a smaller p-value , also suggesting that the residuals are white noise.
Fitting the ARIMA model with the existing data
The residual errors seem fine with near zero mean and uniform variance. Let’s plot the actuals against the fitted values
Convert model and time series to dataframe for plotting
daily_hosp_hb_timeseries_data <- fortify(daily_hosp_hb_timeseries) %>%
clean_names() %>%
remove_rownames %>%
rename (date = index,
hosp = data)%>%
mutate(index = seq(1:nrow(daily_hosp_hb_timeseries)))
arima_fit_resid <- ts(daily_hosp_hb_timeseries[1:nrow(daily_hosp_hb_timeseries)]) - resid(auto_arima_fit_hosp)
arima_fit_data <- fortify(arima_fit_resid) %>%
clean_names() %>%
mutate(data = round(data,2))
fit_existing_data <- daily_hosp_hb_timeseries_data %>%
inner_join(arima_fit_data, by = c("index"))
Plotting the series along with the fitted values
fit_existing_hosp_plot <- fit_existing_data %>%
mutate (date = as.Date(date)) %>%
ggplot()+
aes(x=date, y = hosp)+
geom_line(color ="#5ab4ac")+
geom_line(aes(x= date, y = data), colour = "red" )+
xlab("Month") +
ylab("Patient Hospitalised")+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
ggtitle("Fitting the ARIMA model with existing data") +
scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
#scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))+
color_theme()
ggplotly(fit_existing_hosp_plot)
Step 6 Forecast using the model
Data Preparation :
forecast_model <- forecast(auto_arima_fit_hosp,level = c(80, 95), h = 30)
#Convert the model to dataframe for plotting
forecast_model_data <- fortify(forecast_model) %>%
clean_names() %>%
mutate(data = round(data,2),
fitted= round(fitted,2)) %>%
mutate (lo_80 = ifelse(lo_80 < 0,0,lo_80),
lo_95 = ifelse(lo_95 < 0,0,lo_95)
)
forecast_start_date <- as.Date(max(daily_hosp_hb_timeseries_data$date)+1)
forecast_end_date <- as.Date(forecast_start_date+29)
forecast_data <- forecast_model_data %>%
filter(!(is.na(point_forecast))) %>%
mutate(date = seq(forecast_start_date,forecast_end_date, by =1)) %>%
select(-data,-fitted, -index)
fitted_data <- forecast_model_data %>%
filter(!(is.na(data))) %>%
inner_join(daily_hosp_hb_timeseries_data, by = c("index")) %>%
mutate(date = as.Date(date)) %>%
select(date, data, fitted)
Plotting the Vaccination series plus the forecast and 80 - 95% prediction intervals
annotation <- data.frame(
x = c(as.Date("13-08-2021","%d-%m-%Y"),as.Date("31-10-2021","%d-%m-%Y")),
y = c(180,200),
label = c("PAST", "FUTURE")
)
#Time series plots for the next 60 days according to best ARIMA models with 80%–95% CI.
forecast_data_hosp_plot <-fitted_data %>%
ggplot()+
geom_line(aes(x= date, y = data), color = "#5ab4ac")+
#geom_line(aes(x= date, y = fitted), colour = "red" )+
geom_line(aes(x= date, y =point_forecast), color ="blue", size = 0.5,
data = forecast_data )+
geom_ribbon(aes(x = date, y = point_forecast, ymin = lo_80, ymax = hi_80),
data = forecast_data, alpha = 0.3, fill = "green")+
geom_ribbon(aes(x = date, y = point_forecast, ymin = lo_95, ymax = hi_95),
data = forecast_data, alpha = 0.1)+
geom_vline(aes(xintercept=as.numeric(max(date))),color="#f1a340", linetype="dashed",data = fitted_data)+
ggtitle("Projection of Hospitalisation") +
xlab("Month") +
ylab("Patient Hospitalised")+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
color_theme()+
scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
geom_text(data=annotation,
aes( x=x, y=y, label=label),
color="blue",
size=4 )
ggplotly(forecast_data_hosp_plot)
LS0tDQp0aXRsZTogIkZvcmVjYXN0IG9uIEhvc3BpdGFsaXNhdGlvbiAoQVJJTUEgTW9kZWxsaW5nKSINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KZWRpdG9yX29wdGlvbnM6IA0KICBtYXJrZG93bjogDQogICAgd3JhcDogc2VudGVuY2UNCi0tLQ0KDQoqKkRhdGEgUHJlcGFyYXRpb24qKg0KDQpgYGB7cn0NCiNGb3IgZm9yZWNhc3RpbmcsIHdlIGNob3NlIHRoZSBsYXRlc3QgZGF0YQ0KdHJlbmRfaG9zcF9oYiA8LSB0cmVuZF9oYl9kYWlseSAlPiUgDQogIGZpbHRlciAoaGJfbmFtZSA9PSAiU2NvdGxhbmQiKSAlPiUgDQogIGZpbHRlcihkYXRlID49IjIwMjEtMDYtMDEiKSAlPiUgDQogIGZpbHRlcighKGlzLm5hKGhvc3BpdGFsX2FkbWlzc2lvbnMpKSkgJT4lIA0KICBzZWxlY3QoZGF0ZSwgaG9zcGl0YWxfYWRtaXNzaW9ucykNCg0KIyBDb252ZXJ0IGl0IGludG8gYSB0aW1lIHNlcmllcw0KZGFpbHlfaG9zcF9oYl96b28gPC0gem9vKHRyZW5kX2hvc3BfaGIkaG9zcGl0YWxfYWRtaXNzaW9ucywgDQogICAgICAgICAgIG9yZGVyLmJ5PWFzLkRhdGUodHJlbmRfaG9zcF9oYiRkYXRlLCBmb3JtYXQ9JyVtLyVkLyVZJykpDQoNCiMgQ29udmVydCBpdCBpbnRvIGEgdGltZSBzZXJpZXMNCmRhaWx5X2hvc3BfaGJfdGltZXNlcmllcyA8LSAgdGltZVNlcmllczo6YXMudGltZVNlcmllcyhkYWlseV9ob3NwX2hiX3pvbykNCmBgYA0KDQojIyBTdGVwIDEgOiBWaXN1YWxpemUgdGhlIHRpbWUgc2VyaWVzDQoNCmBgYHtyfQ0Kb3JpZ2luYWxfc2VyaWVzPC1hdXRvcGxvdChkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMsIHRzLmNvbG91ciA9ICcjNWFiNGFjJykrDQogIHhsYWIoIk1vbnRoIikgKyANCiAgeWxhYigiTnVtYmVyIG9mIFBlb3BsZSBob3NwaXRhbGlzZWQiKSsNCiAgI3NjYWxlX3hfZGF0ZShicmVha3MgPSAiMSBtb250aCIsIGRhdGVfbGFiZWxzID0gIiViIC0gJXkiICkrDQogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gNDUsIHZqdXN0ID0gMSwgaGp1c3Q9MSkpKw0KICBnZ3RpdGxlKCJUcmVuZCBvbiBIb3NwaXRhbGlzYXRpb24iKSArDQogIGNvbG9yX3RoZW1lKCkNCg0KZ2dwbG90bHkob3JpZ2luYWxfc2VyaWVzKQ0KYGBgDQoNCiMjIFN0ZXAgMiA6IElkZW50aWZpY2F0aW9uIG9mIG1vZGVsIDogKEZpbmRpbmcgZDopDQoNCklkZW50aWZ5IHdoZXRoZXIgdGhlIHRpbWUgc2VyaWVzIGlzIHN0YXRpb25hcnkgLyBub24gc3RhdGlvbmFyeSB3ZSBjYW4gdXNlIEFERiBBdWdtZW50ZWQgRGlja2V5LUZ1bGxlciB0ZXN0DQoNCmBgYHtyfQ0KYWRmX3Rlc3RfaG9zcCA8LSBhZGYudGVzdChkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMpDQphZGZfdGVzdF9ob3NwDQpgYGANCg0KVGhlIHRpbWUgc2VyaWVzIGlzIG5vdCBzdGF0aW9uYXJ5IHNpbmNlIHdlIGhhdmUgYSBoaWdoIHAtdmFsdWUgKHAtdmFsdWUgbXVzdCBiZSBcPCAwLjA1KS4NClNvIHdlIGFwcGx5IGRpZmZlcmVuY2UNCg0KYGBge3J9DQpmaXJzdF9kaWZmX2hvc3A8LSBkaWZmKGRhaWx5X2hvc3BfaGJfdGltZXNlcmllcykNCmFkZl90ZXN0MV9ob3NwIDwtIGFkZi50ZXN0KG5hLm9taXQoZmlyc3RfZGlmZl9ob3NwKSkNCmFkZl90ZXN0MV9ob3NwDQpgYGANCg0KQ3JlYXRlIGEgZGF0YWZyYW1lIHRvIGNvbXBhcmUNCg0KYGBge3J9DQphZGZfZGF0YV9ob3NwIDwtIGRhdGEuZnJhbWUoRGF0YSA9IGMoIk9yaWdpbmFsIiwgIkZpcnN0LU9yZGVyZWQiKSwNCiAgICAgICAgICAgICAgICAgICAgICAgRGlja2V5X0Z1bGxlciA9IGMoYWRmX3Rlc3RfaG9zcCRzdGF0aXN0aWMsIGFkZl90ZXN0MV9ob3NwJHN0YXRpc3RpYyksDQogICAgICAgICAgICAgICAgICAgICAgIHBfdmFsdWUgPSBjKGFkZl90ZXN0X2hvc3AkcC52YWx1ZSxhZGZfdGVzdDFfaG9zcCRwLnZhbHVlKSkNCmFkZl9kYXRhX2hvc3ANCmBgYA0KDQpJbml0aWFsbHkgdGhlIHAtdmFsdWUgaXMgaGlnaCB3aGljaCBpbmRpY2F0ZXMgdGhhdCB0aGUgVGltZSBTZXJpZXMgaXMgbm90IHN0YXRpb25hcnkuDQpTbyB3ZSBhcHBseSBkaWZmZXJlbmNlIDEgdGltZS4NCkFmdGVyIHRoZSBmaXJzdCBkaWZmZXJlbmNlLCB0aGUgcC12YWx1ZSBcPCBzaWduaWZpY2FuY2UgbGV2ZWwgKDAuMDUpIFNvIHdlIGNhbiBjb25jbHVkZSB0aGF0IHRoZSBkaWZmZXJlbmNlIGRhdGEgYXJlIHN0YXRpb25hcnkuDQpTbyBkaWZmZXJlbmNlIChkID0gMSkNCg0KT3RoZXIgbWV0aG9kOg0KDQpgYGB7cn0NCm5kaWZmcyhkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMpDQpgYGANCg0KTGV0J3MgcGxvdCB0aGUgRmlyc3QgT3JkZXIgRGlmZmVyZW5jZSBTZXJpZXMNCg0KT3JkZXIgb2YgZmlyc3QgZGlmZmVyZW5jZQ0KDQpgYGB7cn0NCnA8LSBhdXRvcGxvdChmaXJzdF9kaWZmX2hvc3AsIHRzLmNvbG91ciA9ICcjNWFiNGFjJykgKw0KICB4bGFiKCJNb250aCIpICsgDQogIHlsYWIoIkhPU1BJVEFMSVpBVElPTiIpKw0KICMgc2NhbGVfeF9kYXRlKGJyZWFrcyA9ICIxIG1vbnRoIiwgZGF0ZV9sYWJlbHMgPSAiJWIgLSAleSIgKSsNCiAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA0NSwgdmp1c3QgPSAxLCBoanVzdD0xKSkrDQogIGdndGl0bGUoIkZpcnN0LU9yZGVyIERpZmZlcmVuY2UgU2VyaWVzIikgKw0KICBjb2xvcl90aGVtZSgpDQoNCmdncGxvdGx5KHApDQpgYGANCg0KIyMgU3RlcCAzOiBFc3RpbWF0ZSB0aGUgcGFyYW1ldGVycyAoRmluZGluZyBwIGFuZCBxKQ0KDQpGb3Igb3VyIG1vZGVsIEFSSU1BIChwLGQscSksIHdlIGZvdW5kIGQgPSAxLCB0aGUgbmV4dCBzdGVwIGlzIHRvIGdldCB0aGUgdmFsdWVzIG9mIHAgYW5kIHEsIHRoZSBvcmRlciBvZiBBUiBhbmQgTUEgcGFydC4NClBsb3QgQUNGIGFuZCBQQUNGIGNoYXJ0cyB0byBpZGVudGlmeSBxIGFuZCBwIHJlc3BlY3RpdmVseS4NCg0KYGBge3J9DQpwYXIobWZyb3c9YygyLDIpKQ0KYWNmX2hvc3AgIDwtIGFjZjEoZGFpbHlfaG9zcF9oYl90aW1lc2VyaWVzLCBjb2w9Mjo3LCBsd2Q9NCkNCnBhY2ZfaG9zcCA8LSBhY2YxKGRhaWx5X2hvc3BfaGJfdGltZXNlcmllcywgIHBhY2YgPSBUUlVFLCBjb2w9Mjo3LCBsd2Q9NCkNCmFjZl9ob3NwICA8LSBhY2YxKGZpcnN0X2RpZmZfaG9zcCwgY29sPTI6NywgbHdkPTQpDQpwYWNmX2hvc3AgPC0gYWNmMShmaXJzdF9kaWZmX2hvc3AsICBwYWNmID0gVFJVRSwgY29sPTI6NywgbHdkPTQpDQpgYGANCg0KVGhlIEFDRiBhbmQgUEFDRiBwbG90cyBvZiB0aGUgZGlmZmVyZW5jZWQgZGF0YSBzaG93IHRoZSBmb2xsb3dpbmcgcGF0dGVybnM6DQoNClRoZSBBQ0YgaXMgc2ludXNvaWRhbCBhbmQgdGhlcmUgaXMgYSBzaWduaWZpY2FudCBzcGlrZSBhdCBsYWcgMyBpbiB0aGUgUEFDRiBTbyB0aGUgZGF0YSBtYXkgZm9sbG93IGFuIEFSKDMpIG1vZGVsDQoNClRoZSBQQUNGIGlzIHNpbnVzb2lkYWwgYW5kIHRoZXJlIGlzIGEgc2lnbmlmaWNhbnQgc3Bpa2UgYXQgbGFnIDIgaW4gdGhlIEFDRiBTbyB0aGUgZGF0YSBtYXkgZm9sbG93IGFuIE1BKDIpIG1vZGVsDQoNClNvIHdlIHByb3Bvc2UgdGhyZWUgQVJNQSBtb2RlbHMgZm9yIHRoZSBkaWZmZXJlbmNlZCBkYXRhOiBBUk1BKHAscSkgQVJNQSgzLDIpLCBBUk1BKDMsMCkgYW5kIEFSTUEoMCwyKS4NCg0KVGhhdCBpcywgZm9yIHRoZSBvcmlnaW5hbCB0aW1lIHNlcmllcywgd2UgcHJvcG9zZSB0aHJlZSBBUklNQSBtb2RlbHMsQVJJTUEocCxkLHEpIEFSSU1BKDMsMSwyKSwgQVJJTUEoMywxLDApIGFuZCBBUk1BKDAsMSwyKS4NCg0KIyMgU3RlcCA0OiBCdWlsZCB0aGUgQVJJTUEgbW9kZWwNCg0KIyMjICoqTWFudWFsIEFSSU1BOioqDQoNCmBgYHtyfQ0KYXJpbWFfZml0X2hvc3BfMSA9IEFyaW1hKGRhaWx5X2hvc3BfaGJfdGltZXNlcmllcywgb3JkZXIgPSBjKDMsMSwyKSkNCmFyaW1hX2ZpdF9ob3NwXzIgPSBBcmltYShkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMsIG9yZGVyID0gYygzLDEsMCkpDQphcmltYV9maXRfaG9zcF8zID0gQXJpbWEoZGFpbHlfaG9zcF9oYl90aW1lc2VyaWVzLCBvcmRlciA9IGMoMCwxLDIpKQ0KYGBgDQoNCmBgYHtyfQ0Kc3VtbWFyeShhcmltYV9maXRfaG9zcF8xKQ0Kc3VtbWFyeShhcmltYV9maXRfaG9zcF8yKQ0Kc3VtbWFyeShhcmltYV9maXRfaG9zcF8zKQ0KYGBgDQoNCmBgYHtyfQ0KDQp0ZXhyZWc6OnNjcmVlbnJlZyhsaXN0KGFyaW1hX2ZpdF9ob3NwXzEsIGFyaW1hX2ZpdF9ob3NwXzIsIGFyaW1hX2ZpdF9ob3NwXzMpLA0KICAgICAgICAgICAgICAgIGN1c3RvbS5tb2RlbC5uYW1lcyA9YygiQVJJTUEoMywxLDIpIiwiQVJJTUEoMywxLDApIiwiQVJJTUEoMCwxLDQpIiksDQogICAgICAgICAgICAgICAgY2VudGVyID0gVFJVRSwNCiAgICAgICAgICAgICAgICB0YWJsZSA9IEZBTFNFKQ0KYGBgDQoNCmBgYHtyfQ0KI0Z1bmN0aW9uIGZvciBBdXRvbWF0ZWQgQVJJTUENCg0KDQphdXRvX2FyaW1hX2ZpdF9ob3NwIDwtIGF1dG8uYXJpbWEobGFnKGRhaWx5X2hvc3BfaGJfdGltZXNlcmllcyxrID1sYWdfdmFsdWUgKSwNCiAgICAgICAgICAgICAgICAgIHNlYXNvbmFsPUZBTFNFLA0KICAgICAgICAgICAgICAgICAgc3RlcHdpc2U9RkFMU0UsDQogICAgICAgICAgICAgICAgICBhcHByb3hpbWF0aW9uPUZBTFNFLA0KICAgICAgICAgICAgICAgICAgdHJhY2UgPSBUUlVFDQogICAgICAgICAgICAgICAgICApDQpgYGANCg0KIyMjIEF1dG9tYXRlZCBBUklNQQ0KDQpgYGB7cn0NCiNMYWcgaXMgdXNlZCB0byBiZXN0IGZpdCB0aGUgbW9kZWwNCmF1dG9fYXJpbWFfZml0X2hvc3AgPC0gYXV0by5hcmltYShsYWcoZGFpbHlfaG9zcF9oYl90aW1lc2VyaWVzKSwNCiAgICAgICAgICAgICAgICAgIHNlYXNvbmFsPUZBTFNFLA0KICAgICAgICAgICAgICAgICAgc3RlcHdpc2U9RkFMU0UsDQogICAgICAgICAgICAgICAgICBhcHByb3hpbWF0aW9uPUZBTFNFLA0KICAgICAgICAgICAgICAgICAgdHJhY2UgPSBUUlVFDQogICAgICAgICAgICAgICAgICApDQphdXRvX2FyaW1hX2ZpdF9ob3NwDQpgYGANCg0KQXV0b21hdGVkIEFSSU1BIGNvbmZpcm1zIHRoYXQgdGhlIEFSSU1BKDMsMSwyKSBzZWVtcyBnb29kIGJhc2VkIG9uIEFJQw0KDQpgYGB7cn0NCmNvZWY8LWxtdGVzdDo6Y29lZnRlc3QoYXV0b19hcmltYV9maXRfaG9zcCkNCmNvZWYNCmBgYA0KDQpBbGwgY29lZmZpY2llbnRzIGFyZSBzaWduaWZpY2FudCBleGNlcHQgYXIzLg0KDQojIyMgTW9kZWwgU2VsZWN0aW9uIENyaXRlcmlhIDoNCg0KQVJJTUEgbW9kZWxzIHdpdGggbWluaW11bSBBSUMsIFJNU0UgYW5kIE1BUEUgY3JpdGVyaWEgd2VyZSBjaG9zZW4gYXMgdGhlIGJlc3QgbW9kZWxzLg0KQmFzZWQgb24gQWthaWtlIEluZm9ybWF0aW9uIENyaXRlcmlvbiAoQUlDKSBhYm92ZSwgYW4gQVJJTUEoMywgMSwgMikgbW9kZWwgc2VlbXMgYmVzdC4NCg0KIyMgU3RlcCA1OiBDaGVjayBmb3IgRGlhZ25vc3RpY3MNCg0KTGV0J3MgcGxvdCB0aGUgZGlhZ25vc3RpY3Mgd2l0aCB0aGUgcmVzdWx0cyB0byBtYWtlIHN1cmUgdGhlIG5vcm1hbGl0eSBhbmQgY29ycmVsYXRpb24gYXNzdW1wdGlvbnMgZm9yIHRoZSBtb2RlbCBob2xkLg0KSWYgdGhlIHJlc2lkdWFscyBsb29rIGxpa2Ugd2hpdGUgbm9pc2UsIHByb2NlZWQgd2l0aCBmb3JlY2FzdCBhbmQgcHJlZGljdGlvbiwgb3RoZXJ3aXNlIHJlcGVhdCB0aGUgbW9kZWwgYnVpbGRpbmcuDQoNCmBgYHtyfQ0KcmVzIDwtY2hlY2tyZXNpZHVhbHMoYXV0b19hcmltYV9maXRfaG9zcCwgdGhlbWUgPSBjb2xvcl90aGVtZSgpKQ0KcmVzDQpgYGANCg0KVGhlIEFDRiBwbG90IG9mIHRoZSByZXNpZHVhbHMgZnJvbSB0aGUgQVJJTUEoMywxLDIpIG1vZGVsIHNob3dzIHRoYXQgYWxsIGF1dG8gY29ycmVsYXRpb25zIGFyZSBhbG1vc3Qgd2l0aGluIHRoZSB0aHJlc2hvbGQgbGltaXRzLCB3aXRoIHJlc2lkdWFscy4NCkEgcG9ydG1hbnRlYXUgdGVzdCAoTGp1bmctQm94IHRlc3QpIHJldHVybnMgYSBzbWFsbGVyIHAtdmFsdWUgLCBhbHNvIHN1Z2dlc3RpbmcgdGhhdCB0aGUgcmVzaWR1YWxzIGFyZSB3aGl0ZSBub2lzZS4NCg0KRml0dGluZyB0aGUgQVJJTUEgbW9kZWwgd2l0aCB0aGUgZXhpc3RpbmcgZGF0YQ0KDQpUaGUgcmVzaWR1YWwgZXJyb3JzIHNlZW0gZmluZSB3aXRoIG5lYXIgemVybyBtZWFuIGFuZCB1bmlmb3JtIHZhcmlhbmNlLg0KTGV0J3MgcGxvdCB0aGUgYWN0dWFscyBhZ2FpbnN0IHRoZSBmaXR0ZWQgdmFsdWVzDQoNCioqQ29udmVydCBtb2RlbCBhbmQgdGltZSBzZXJpZXMgdG8gZGF0YWZyYW1lIGZvciBwbG90dGluZyoqDQoNCmBgYHtyfQ0KZGFpbHlfaG9zcF9oYl90aW1lc2VyaWVzX2RhdGEgPC0gZm9ydGlmeShkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMpICU+JSANCiAgY2xlYW5fbmFtZXMoKSAlPiUgDQogIHJlbW92ZV9yb3duYW1lcyAlPiUgDQogIHJlbmFtZSAoZGF0ZSA9IGluZGV4LA0KICAgICAgICAgIGhvc3AgPSBkYXRhKSU+JSANCiAgbXV0YXRlKGluZGV4ID0gc2VxKDE6bnJvdyhkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMpKSkNCiAgDQphcmltYV9maXRfcmVzaWQgPC0gdHMoZGFpbHlfaG9zcF9oYl90aW1lc2VyaWVzWzE6bnJvdyhkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXMpXSkgLSByZXNpZChhdXRvX2FyaW1hX2ZpdF9ob3NwKQ0KDQphcmltYV9maXRfZGF0YSA8LSBmb3J0aWZ5KGFyaW1hX2ZpdF9yZXNpZCkgJT4lIA0KICBjbGVhbl9uYW1lcygpICU+JSANCiAgbXV0YXRlKGRhdGEgPSByb3VuZChkYXRhLDIpKQ0KDQpmaXRfZXhpc3RpbmdfZGF0YSA8LSBkYWlseV9ob3NwX2hiX3RpbWVzZXJpZXNfZGF0YSAlPiUgDQogIGlubmVyX2pvaW4oYXJpbWFfZml0X2RhdGEsIGJ5ID0gYygiaW5kZXgiKSkNCmBgYA0KDQoqKlBsb3R0aW5nIHRoZSBzZXJpZXMgYWxvbmcgd2l0aCB0aGUgZml0dGVkIHZhbHVlcyoqDQoNCmBgYHtyfQ0KZml0X2V4aXN0aW5nX2hvc3BfcGxvdCA8LSBmaXRfZXhpc3RpbmdfZGF0YSAlPiUgDQogIG11dGF0ZSAoZGF0ZSA9IGFzLkRhdGUoZGF0ZSkpICU+JSANCiAgZ2dwbG90KCkrDQogIGFlcyh4PWRhdGUsIHkgPSBob3NwKSsNCiAgZ2VvbV9saW5lKGNvbG9yID0iIzVhYjRhYyIpKw0KICBnZW9tX2xpbmUoYWVzKHg9IGRhdGUsIHkgPSBkYXRhKSwgY29sb3VyID0gInJlZCIgKSsNCiAgeGxhYigiTW9udGgiKSArIA0KICB5bGFiKCJQYXRpZW50IEhvc3BpdGFsaXNlZCIpKw0KICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDQ1LCB2anVzdCA9IDEsIGhqdXN0PTEpKSsNCiAgZ2d0aXRsZSgiRml0dGluZyB0aGUgQVJJTUEgbW9kZWwgd2l0aCBleGlzdGluZyBkYXRhIikgKw0KICBzY2FsZV94X2RhdGUoYnJlYWtzID0gIjEgbW9udGgiLCBkYXRlX2xhYmVscyA9ICIlYiAtICV5IiApKw0KICAjc2NhbGVfeV9jb250aW51b3VzKGxhYmVscyA9IHNjYWxlczo6dW5pdF9mb3JtYXQodW5pdCA9ICJNIiwgc2NhbGUgPSAxZS02KSkrDQogIGNvbG9yX3RoZW1lKCkNCg0KZ2dwbG90bHkoZml0X2V4aXN0aW5nX2hvc3BfcGxvdCkNCmBgYA0KDQojIyBTdGVwIDYgRm9yZWNhc3QgdXNpbmcgdGhlIG1vZGVsDQoNCioqRGF0YSBQcmVwYXJhdGlvbiA6KioNCg0KYGBge3J9DQpmb3JlY2FzdF9tb2RlbCA8LSBmb3JlY2FzdChhdXRvX2FyaW1hX2ZpdF9ob3NwLGxldmVsID0gYyg4MCwgOTUpLCBoID0gMzApIA0KDQojQ29udmVydCB0aGUgbW9kZWwgdG8gZGF0YWZyYW1lIGZvciBwbG90dGluZw0KDQpmb3JlY2FzdF9tb2RlbF9kYXRhIDwtIGZvcnRpZnkoZm9yZWNhc3RfbW9kZWwpICU+JSANCiAgY2xlYW5fbmFtZXMoKSAlPiUgDQogIG11dGF0ZShkYXRhID0gcm91bmQoZGF0YSwyKSwNCiAgICAgICAgIGZpdHRlZD0gcm91bmQoZml0dGVkLDIpKSAlPiUgDQogIG11dGF0ZSAobG9fODAgPSBpZmVsc2UobG9fODAgPCAwLDAsbG9fODApLA0KICAgICAgICAgIGxvXzk1ID0gaWZlbHNlKGxvXzk1IDwgMCwwLGxvXzk1KQ0KICApDQoNCmZvcmVjYXN0X3N0YXJ0X2RhdGUgPC0gYXMuRGF0ZShtYXgoZGFpbHlfaG9zcF9oYl90aW1lc2VyaWVzX2RhdGEkZGF0ZSkrMSkNCmZvcmVjYXN0X2VuZF9kYXRlIDwtIGFzLkRhdGUoZm9yZWNhc3Rfc3RhcnRfZGF0ZSsyOSkNCg0KZm9yZWNhc3RfZGF0YSA8LSBmb3JlY2FzdF9tb2RlbF9kYXRhICU+JSANCiAgZmlsdGVyKCEoaXMubmEocG9pbnRfZm9yZWNhc3QpKSkgJT4lIA0KICBtdXRhdGUoZGF0ZSA9IHNlcShmb3JlY2FzdF9zdGFydF9kYXRlLGZvcmVjYXN0X2VuZF9kYXRlLCBieSA9MSkpICU+JSANCnNlbGVjdCgtZGF0YSwtZml0dGVkLCAtaW5kZXgpICANCg0KZml0dGVkX2RhdGEgPC0gZm9yZWNhc3RfbW9kZWxfZGF0YSAlPiUgDQogIGZpbHRlcighKGlzLm5hKGRhdGEpKSkgJT4lIA0KICBpbm5lcl9qb2luKGRhaWx5X2hvc3BfaGJfdGltZXNlcmllc19kYXRhLCBieSA9IGMoImluZGV4IikpICU+JSANCiAgbXV0YXRlKGRhdGUgPSBhcy5EYXRlKGRhdGUpKSAlPiUgDQpzZWxlY3QoZGF0ZSwgZGF0YSwgZml0dGVkKSANCg0KYGBgDQoNCioqUGxvdHRpbmcgdGhlIFZhY2NpbmF0aW9uIHNlcmllcyBwbHVzIHRoZSBmb3JlY2FzdCBhbmQgODAgLSA5NSUgcHJlZGljdGlvbiBpbnRlcnZhbHMqKg0KDQpgYGB7cn0NCg0KYW5ub3RhdGlvbiA8LSBkYXRhLmZyYW1lKA0KICAgeCA9IGMoYXMuRGF0ZSgiMTMtMDgtMjAyMSIsIiVkLSVtLSVZIiksYXMuRGF0ZSgiMzEtMTAtMjAyMSIsIiVkLSVtLSVZIikpLA0KICAgeSA9IGMoMTgwLDIwMCksDQogICBsYWJlbCA9IGMoIlBBU1QiLCAiRlVUVVJFIikNCikNCg0KI1RpbWUgc2VyaWVzIHBsb3RzIGZvciB0aGUgbmV4dCA2MCBkYXlzIGFjY29yZGluZyB0byBiZXN0IEFSSU1BIG1vZGVscyB3aXRoIDgwJeKAkzk1JSBDSS4NCmZvcmVjYXN0X2RhdGFfaG9zcF9wbG90IDwtZml0dGVkX2RhdGEgJT4lIA0KICBnZ3Bsb3QoKSsNCiAgZ2VvbV9saW5lKGFlcyh4PSBkYXRlLCB5ID0gZGF0YSksIGNvbG9yID0gIiM1YWI0YWMiKSsNCiAgI2dlb21fbGluZShhZXMoeD0gZGF0ZSwgeSA9IGZpdHRlZCksIGNvbG91ciA9ICJyZWQiICkrDQogIGdlb21fbGluZShhZXMoeD0gZGF0ZSwgeSA9cG9pbnRfZm9yZWNhc3QpLCBjb2xvciA9ImJsdWUiLCBzaXplID0gMC41LA0KICAgICAgICAgICAgIGRhdGEgPSBmb3JlY2FzdF9kYXRhICkrDQogIGdlb21fcmliYm9uKGFlcyh4ID0gZGF0ZSwgeSA9IHBvaW50X2ZvcmVjYXN0LCB5bWluID0gbG9fODAsIHltYXggPSBoaV84MCksIA0KICAgICAgICAgICAgICBkYXRhID0gZm9yZWNhc3RfZGF0YSwgYWxwaGEgPSAwLjMsIGZpbGwgPSAiZ3JlZW4iKSsNCiAgZ2VvbV9yaWJib24oYWVzKHggPSBkYXRlLCB5ID0gcG9pbnRfZm9yZWNhc3QsIHltaW4gPSBsb185NSwgeW1heCA9IGhpXzk1KSwgDQogICAgICAgICAgICAgIGRhdGEgPSBmb3JlY2FzdF9kYXRhLCBhbHBoYSA9IDAuMSkrDQogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQ9YXMubnVtZXJpYyhtYXgoZGF0ZSkpKSxjb2xvcj0iI2YxYTM0MCIsIGxpbmV0eXBlPSJkYXNoZWQiLGRhdGEgPSBmaXR0ZWRfZGF0YSkrDQogIGdndGl0bGUoIlByb2plY3Rpb24gb2YgSG9zcGl0YWxpc2F0aW9uIikgKw0KICB4bGFiKCJNb250aCIpICsgDQogIHlsYWIoIlBhdGllbnQgSG9zcGl0YWxpc2VkIikrDQogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gNDUsIHZqdXN0ID0gMSwgaGp1c3Q9MSkpKw0KICBjb2xvcl90aGVtZSgpKw0KICBzY2FsZV94X2RhdGUoYnJlYWtzID0gIjEgbW9udGgiLCBkYXRlX2xhYmVscyA9ICIlYiAtICV5IiApKw0KICAgZ2VvbV90ZXh0KGRhdGE9YW5ub3RhdGlvbiwgDQogICAgICAgICAgICAgYWVzKCB4PXgsIHk9eSwgbGFiZWw9bGFiZWwpLCAgICAgICAgICAgICAgICAgIA0KICAgICAgICAgICAgY29sb3I9ImJsdWUiLCANCiAgICAgICAgICAgIHNpemU9NCApDQogIA0KICBnZ3Bsb3RseShmb3JlY2FzdF9kYXRhX2hvc3BfcGxvdCkNCiAgDQpgYGANCg==